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Abstract 

Five models for matrix damage in fiber reinforced laminates previously presented 
in the literature are evaluated for matrix dominated loading conditions under plane 
stress (i.e. tension transverse to fibers combined with in-plane shear) and are com- 
pared both qualitatively and quantitatively. The emphasis of this study is on mod- 
eling the response of embedded plies to homogeneous stress states rather than lo- 
calized failure caused by stress concentrations (e.g. at notches). Three of the models 
are specifically designed for modeling the non-linear response due to distributed ma- 
trix cracking under homogeneous loading, and also account for non-linear (shear) 
behavior prior to the onset of cracking. The remaining two models are localized 
damage models intended for predicting local failure at stress concentrations. These 
last two models are not strictly applicable to the test cases considered here but, nev- 
ertheless, are also included in this study for comparison and to point out qualitative 
differences in the various modeling strategies. 

In this study, the modeling approaches of distributed vs. localized cracking as 
well as the different formulations of damage initiation and damage progression are 
discussed. In addition, the following issues are addressed: 

• stress-strain response of an embedded ply under transverse tension, 

• damage initiation, 

• anisotropy of property degradation due to damage, 

• non-linearity prior to damage onset (plasticity). 


Key words: Fiber Reinforced Laminates, Polymer Matrix Composites, 
Computational Mechanics, Non-Linear Modeling, Continuum Damage, Plasticity. 



Notation 


Indices: 

1, 2, 3 ... ply coordinates (1 - fiber, 2 - transverse in-plane, 3 - out-of-plane direction) 

I, n, t ... fracture plane coordinates (I - fiber, n - normal, t - transverse direction) 
x, y, z ... global coordinates (x - loading direction, y - transverse to load, in plane, 
z - out of plane) 


Greek letters: 

P 

lij 

pi 

lij 

A U /width 
5eq 

S° eq 

su 
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£-ii 

£ pl 

-eq 

Vl 
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P21 

er 
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'U 
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_eff 
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eff 
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eff 


°fp 
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lay-up angle for off-axis plies 

engineering shear strain on plane i in direction j 

plastic shear strain on plane i in direction j 

energy dissipated by the creation of a new crack 

equivalent displacement (ABAQUS model) 

equivalent displacement at damage onset (ABAQUS model) 

equivalent displacement at final failure (ABAQUS model) 

strain tensor 

normal strain component on plane i in direction i 
plastic strain tensor 

equivalent strain ( constant stress model) 

slope parameter for LaRC 04 criterion 

damage evolution parameter {Mori- Tanaka model) 

major in-plane Poisson’s ratio 

minor in-plane Poisson’s ratio (zz 2 i = zq 2 Ei/E\) 

ply stress tensor 

stress component on plane i in direction j 

stress component at damage onset on plane i in direction j 

effective stress tensor (ABAQUS model) 

effective stress tensor ( constant stress and periodic crack model) 

magnitude of effective stress tensor 

effective stress component on plane i in direction j 

fracture plane stress (magnitude of fracture plane traction vector) 

fracture plane stress at damage onset 

fracture plane angle 

shear response factor ( periodic crack model) 
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Roman letters: 


CD 

CD sat 

d { 

dm 

d s 

E d 

■^Lam 

Ei 

Ef 

El 

jfjsec 

fl2 

h 

FI 

C\ 2 

mo 

Lt 12 

G sec 
12 

dr Ir.ply 
G lie, ply 
^c,ply 
/iges 

hyo 

k 

L 

Lc\i 

^sat 


n 

Pn 

S 

S 1F 

u dis 

jj ~ sat 

t 

yt 

J linear 

yt 

J non-linear 


crack density 

crack density at saturation 

damage variable for fiber failure (ABAQUS model) 

damage variable for matrix failure (ABAQUS model) 

damage variable for shear failure (ABAQUS model) 

elasticity tensor of a damaged ply 

elasticity tensor of a laminate ( periodic crack model) 

Young’s modulus of a ply in i-direction 
Young’s modulus of a damaged ply in z-direction 
Young’s modulus of an undamaged ply in transverse direction 
secant Young’s modulus in transverse direction 
inclusion aspect ratio (Mori- Tanaka damage model) 
non-linear function, in-plane shear (constant stress model) 
non-linear function, transverse tension (constant stress model) 
failure index 

in-plane shear modulus of a ply 
shear modulus of an undamaged ply 
secant shear modulus 

critical energy release rate for intra-ply cracking in mode I 

critical energy release rate for intra-ply cracking in mode II 

critical energy release rate for intra-ply cracking (mixed mode) 

half thickness of laminate (periodic crack model) 

half thickness of 90° ply-stack (periodic crack model) 

parameter of shear plasticity law (Mori-Tanaka model) 

half distance between two cracks (periodic crack model) 

characteristic length (exponential softening and ABAQUS model) 

half distance between two cracks at crack saturation 

exponent of shear plasticity law (Mori-Tanaka model) 

slope parameter for Puck failure criterion 

in-situ shear strength 

interlaminar shear strength 

dissipated energy density 

dissipated energy density at crack saturation 

thickness of a cluster of equally oriented plies 

in-situ transverse tensile strength (linear material) 

in-situ transverse tensile strength (non-linear material) 
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1 INTRODUCTION 


Fiber reinforced composites are being increasingly used in applications that 
require lightweight design. In order to achieve further weight reductions with- 
out compromising the reliability of composite parts, accurate prediction of the 
material response including non-linear behavior is beneficial. Models for pre- 
dicting the non-linear response of laminates are typically formulated on the 
ply level and distinguish between various failure modes in order to capture 
the different effects of those failure modes appropriately. 

In polymer matrix composites, the first failure modes that are normally ob- 
served are dominated by the matrix constituent. These failure modes affect 
the response of a ply, and may lead to a non-linear stress-strain relation, but 
in most cases do not lead to complete laminate failure. Depending on the kind 
of loading, the failure modes that can be observed are matrix cracking and 
matrix plasticity. Many models for predicting the non-linear response of an 
embedded ply have been proposed in the literature. While most of the models 
focus on stiffness degradation due to matrix cracking (e.g. 1-9) some also 
include non-linearity caused by matrix plasticity (e.g. 10 - 12). 

In the present work, five damage models that have recently been presented 
in the literature are evaluated [5, 7, 9, 11, 12]. The models were chosen to 
provide a sample of different approaches that are representative of the state 
of the art in damage modeling. The first three damage models include for- 
mulations for non-linearity due to matrix plasticity and to progressive matrix 
cracking [5, 11, 12]. The remaining two models include non-linearity due to 
localized progressive matrix cracking only [7, 9]. A brief description of each 
model is presented, followed by a general discussion on modeling issues, such 
as distributed vs. localized cracking and different formulations of damage ini- 
tiation and damage progression. Finally, predictions from these models using 
the same set of input data are compared, qualitatively and quantitatively, to 
experimental data from the literature, evaluating their ability to accurately 
predict the change in laminate axial stiffness and Poisson’s ratio as a function 
of applied strain. The objective of this comparison is to point out the respec- 
tive features and differences of the various models regarding plasticity, damage 
initiation, and property degradation due to damage, and also to provide some 
guidelines with respect to the applicability of different modeling strategies to 
a given problem. 

A state of plane stress is assumed throughout this study even though some 
of the models can also be applied to triaxial stress states. Furthermore, the 
discussions here are focused on loading scenarios that lead to a combination 
of transverse tension and in-plane shear stresses in each ply. According to ply 
failure criteria that are based on Mohr’s fracture hypothesis for brittle mate- 
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rials such as the Puck criterion [13] and the LaRC04 criterion [14], such stress 
states cause matrix cracks that are perpendicular to the laminate plane. Dam- 
age due to inclined matrix cracks, which occur under transverse compression, 
are not investigated here. Finally, it is noted that the laminates studied here 
were designed to emphasize the effect of matrix damage such that it could 
be studied properly. Laminates for practical applications usually have more 
evenly distributed fiber directions (e.g. quasi-isotropic laminates) and thinner 
ply groups, so the effect of matrix cracking on laminate stiffness is typically 
smaller. 


2 OVERVIEW OF MODELS 


In this section, the damage models considered in this study are summarized. 
The general idea and basic assumptions of the models are explained, rather 
than providing complete model definitions. For detailed information on the 
model formulations, see the respective references. 


2. 1 Damage / plasticity model using Mori- Tanaka approach 


The ‘MoriT’ model proposed by Schuecker et al. [10, 11] is a constitutive 
model specifically for matrix-dominated non-linearities in plane stress. It in- 
cludes non-linearity due to progressive matrix cracking and matrix plasticity, 
while damage due to fiber failure modes is not taken into account. Instead, 
fiber failure is considered as the ultimate failure. The damage and plasticity 
laws used in the model are both formulated with respect to a fracture plane 
that is parallel to the fiber direction and defined by the angle 9f p (see Fig.l). 
The fracture plane orientation is predicted by the Puck failure criterion for 
plane stress states [13, 15, 16]. 

Plastic strain is assumed to be driven by the shear stress in the fracture plane 
(i.e. the projection of the traction vector onto the fracture plane) and to follow 
the direction of the fracture plane’s shear stress vector. Therefore, the only 



Fig. 1. Definition of fracture plane and corresponding coordinate system, l-n-t, with 
regard to the ply coordinate system, 1-2-3, by fracture plane angle, 9{ p . 
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non-zero contributions to plastic strain, when referenced to the fracture plane 
coordinates l-n-t, are the shear strain components 

7ln' = / (crin = 0-12 COS 9 [p ) ^ 

7nt = / (o'nt = -Oil COS 9 ip Sm 9{ p ) . 

For the load cases considered in this study, the Puck failure criterion predicts 
a perpendicular fracture plane (9{ p — 0). This means that the fracture plane 
coordinate system coincides with the ply coordinate system (cf. Fig.l), leaving 
the in-plane shear strain as the only non-zero component of the plastic strain 
tensor ( 7 ^ = y^, 1 ). The relation yjv, = f(ou) is determined by analytical curve 
fitting to experimental data. 

Damage onset and damage evolution are controlled by the failure index com- 
puted from the Puck failure criterion. The stresses used to evaluate the Puck 
criterion are the nominal (homogenized) stresses in the ply and are computed 
as 

<r = E d : (e - e pl ) , (2) 

where e is the given strain tensor, s pl is the plastic strain tensor determined 
by the plasticity law, and E d is the elasticity tensor of the damaged material. 
The elasticity tensor E d is computed by the Mori-Tanaka method [17] using 
penny shaped inclusions which are aligned with the predicted fracture plane. 
The increase of damage is controlled by increasing the volume fraction of voids 
in the Mori-Tanaka formulation according to an empirical evolution law that 
relates the volume fraction to the failure index, FI, of the Puck criterion. Due 
to the formulation of the evolution law proposed in [11], damage onset occurs 
before the failure envelope is reached at the failure index FI — depending 
on the evolution parameter, k. Since FI is a function of the stress state, the 
elasticity tensor (as well as the plastic strain tensor) implicitly depends on 
the stress tensor a. The constitutive equation, Eq. (2), therefore needs to be 
solved iteratively to compute er, e pl (er) and E A (cr) for a given strain state, e. 


2.2 Damage / plasticity model assuming constant stress 


The ‘ConStress ’ model proposed by Pinho et al. [12] is applicable to triaxial 
stress states and includes effects of all currently known failure processes such 
as matrix cracking, matrix plasticity, fiber rupture, fiber kinking, and delami- 
nation. Additionally, it takes into account the influence of hydrostatic stresses 
on elastic properties and fiber kinking. In the present study, only matrix crack- 
ing and matrix plasticity are considered. Furthermore, only plane stress states 
are investigated and the change of elastic properties due to hydrostatic stress 
is therefore neglected. 
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Non-linearity prior to damage onset, which is caused by matrix plasticity, is 
described by the change of secant shear modulus and secant transverse Young’s 
modulus as functions of an equivalent strain (defined here for plane stress) 

ET/E* = / 2 ( £eq ), GT 2 c /G° 12 = / 12 (e eq ) 

/ V^/ 

where e eq = y (e 22 - £ 33 ) 2 + qf 2 . 

The functions / 2 (e e q ) and /i 2 (e e q) need to be determined from experimental 
data of the corresponding uniaxial tests either by point wise interpolation or 
analytical curve fits. 

The failure criterion used to predict the onset of matrix cracking is a modi- 
fied version of the Mohr-Coulomb failure criterion adapted for uni-directional 
composites. Similarly to the Puck failure criterion, this criterion predicts the 
fracture plane orientation, 0f p . in addition to the failure index. Onset and 
propagation of damage are controlled by an effective stress tensor defined as 

<r eff = E° : (e - e pl ) , (4) 

where E° denotes the initial elasticity tensor of the undamaged ply material 
and £ pl is the strain tensor resulting from the plasticity formulation in Eq. (3). 
For perpendicular fracture planes (6{ p — 0), the transverse Young’s modulus 
and the shear modulus are degraded in such a way that the nominal stress 
on the fracture plane, <7f p (defined as the magnitude of the fracture plane’s 
traction vector), remains constant after damage onset, i.e. 

E A 2 = (l-d)El ec , G d l2 = (1 - d)G% 

( j9 

where d — 1 ^ (5) 

<7 ett 

and [, = VV 22) 2 + ( (J i 2 ) 2 ) ^ eff = a/ (Y 22) 2 + (^lf) 2 , 

where < 7 ° denotes the magnitude of the fracture plane traction vector at dam- 
age onset. No additional non-linearity due to plasticity is assumed after dam- 
age onset. Therefore, £ pl does not change during damage accumulation and 
since E° is constant, the effective stress and the ensuing stiffness degradation 
can be computed explicitly from Eqs. (4) and (5) for a given strain state, £. 


2.3 Periodic damage model 


The periodic damage model, ‘PerCrack’ [6, 18], focuses on the prediction 
of transverse matrix cracks. The model is based on an analytical solution 
proposed by Tan and Nuismer [19, 20], which relates the elasticity tensor of 
a laminate consisting of a 90°-ply embedded between two orthotropic sub- 
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laminates, S, to the density of periodically arranged matrix cracks in the 
90°-ply (see Fig. 2). Strictly speaking, the model is therefore restricted to 
laminates that have a (S'/90 n ) s lay-up and in its original version can only 
be applied for either uniaxial tension or in-plane simple shear. The version 
proposed by Camanho and Mayugo [6, 18] combines the Tan and Nuismer 
solution with the LaRC 04 failure criteria for matrix cracking [14] which are 
used to predict damage activation and to control damage evolution under 
multi-axial stress states. The crack density is related to the strain tensor by 
equating the change of elastic strain energy due to stiffness degradation and 
the energy dissipated by the creation of a new crack 

AU/ width = 2/i g esT e:{E Lam (L) - E Lam (L / 2)) : e = 2/), 90 G c , P iy , (6) 

where E La m is the elasticity tensor of the whole laminate as a function of the 
crack spacing, 2 L, which is related to crack density, CD , by 2 L — 1/CD. The 
intra-laminar fracture toughness for a given mode mix, G CiP i y , is determined 
using the mode-mix criterion proposed by Hahn [21]. The only additional 
assumption required to obtain the relation between crack density and the 
strain state is that the multi-axial strain ratio, 7 X y/^ xx (^-direction defined 
as perpendicular to the 90°-ply), remains constant throughout the loading 
history. The LaRC 04 damage activation functions are computed explicitly by 
using the same effective stress tensor as defined in Eq. (4) but with e pl = 0 
since plasticity is not included in this model, i.e. 

o- eff = E° : e . (7) 


To account for the non-linear shear response prior to damage onset, the model 
uses a non-linear elastic relation between in-plane shear stresses and strains 
(without accounting for any influence of transverse stresses) proposed in [22] 

712 = ^ + 7 i N 2 L = ^ + X<t 3 i2 , (8) 

(_ T 12 <-*12 

where 7 is the shear response factor, a material parameter that is determined 
from experimental data. 
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Fig. 2. Periodic cracks in a 90°-ply embedded between orthotropic sub laminates 
[6, 18], 



2-4 Exponential softening model (localized failure) 


In contrast to the previous models, the ‘Exp Soft' model proposed by Maim! 
et al. [7, 8] is designed to model localized failure, i.e. the development of 
a localized crack caused by a stress concentration. It is implemented as a 
user material for the Finite Element Method (FEM) code ABAQUS (Dassault 
Systemes Simulia Corp., Providence, RI, USA) and can be used with plane 
stress (shell) elements and 3D continuum elements. Note however, that only 
the in-plane stress components are used to assess damage evolution, even if 3D 
elements are used. The model includes damage due to tensile and compressive 
failure of both the fibers and the matrix. Here, only degradation due to tensile 
matrix failure is discussed. 


The LaRC 04 failure criteria for transverse matrix cracking [14] are used for 
predicting damage onset and for controlling the progression of damage based 
on the effective stress defined in Eq. (7). After damage onset, the shear modu- 
lus and transverse Young’s modulus are gradually degraded to zero following 
an exponential softening laws (Fig. 3a). In order to produce results that are 
independent of element size, the energy dissipated inside an element needs 
to equal the amount of energy necessary to create exactly one crack through 
that element. This is achieved by introducing a characteristic length of the 
element, L c h, following the crack band model by Bazant [23] and computing 
the energy dissipated per unit volume as 


rrdis 

U 22 

lydis 

U 12 


J a-22d£22 — Gic iP l y /L c h 

J <712^712 = Gllc,ply/^ch 


(9) 


where Gi CjP i y and Gn c>p i y are the intra-laminar fracture toughnesses for mode 
I and mode II, respectively. The softening curves for shear and transverse 
tension are then adjusted such that Eqs. (9) hold. Depending on the chosen 
element size, it is possible that the energy that needs to be dissipated in the 
element exceeds the elastic energy at damage onset, i.e. Gi C)P i y /T c h < \ Y t e% 2 
or Gnc iP i y /T c h < ^ which would lead to “snap-back” in the softening 

curve of the stress-strain law. Therefore, if the element size is chosen too 
large, the strengths Y t and S are reduced appropriately in order to avoid this 
snap-back. 


2.5 ABAQUS damage model for composites (localized failure) 

The ‘ABQ 1 damage model for fiber reinforced composites available in the 
commercial FEM code ABAQUS [9, 24] is in many ways similar to the model 
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described in Section 2.4. It is also conceived as a localizing damage model that 
can handle progressive matrix and fiber failure, each modeled differently for 
tension and compression, and uses Bazant’s crack band model [23] to reduce 
mesh dependency. 

In contrast to the exponential softening model, however, the ABAQUS model 
uses a linear softening law (Fig. 3b). Furthermore, damage onset is predicted 
based on the Hashin failure criterion [25] which is applied to an effective stress 
defined as 

(10) 

where df, d m , and d s are internal variables controlling the degradation of the 
in-plane moduli, 

E* = (l-d t )E lt E d = (1 - d m )E 2 , G d u = (1 - d s )G 12 , . (11) 

For damage due to tensile matrix failure df = 0 and d s — d m , i.e. E d and Gf 2 
follow the same degradation law. Degradation is controlled by an equivalent 
displacement, d eq , which takes the characteristic length into account such that 
the softening curve from damage onset to ultimate failure is linear (Fig. 3b) 
and the energy dissipated by the damage process corresponds to the energy 
required to create one crack inside an element 

U dis L ch = ±a 0 e / eci = G lc , ply , (12) 

where the equivalent stress at damage onset, <Tg q , is a function of the in-plane 
shear and transverse normal stresses and strains. Snap-back is not allowed 
in this model, i.e. 5f q > 5g q always holds true. However, in this model the 
strengths are not adjusted. If the chosen element size is too large, the energy 
dissipated in the element is greater than Gi CjP i y /L c h. 


er = 


l 

1— (if 

o 

0 


0 0 


1 — 0 


0 

1 

1 -rl 




Fig. 3. Softening laws for localized damage models using Bazant’s crack band model; 
(a) exponential softening model [7, 8]; (b) ABAQUS damage model [9, 24]. 
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Table 1 

Material constants of the glass fiber/epoxy material [27, 28, 29] used in the analyses. 


Elastic properties toughnesses slope parameter 


Fi 

e 2 

G u 

V \ 2 

Crlc,ply 

G He, ply 

Pl2 = 

[GPa] 

[GPa] 

[GPa] 


[kJ/m 2 ] 

[kJ/m 2 ] 


44.73 

12.8 

5.8 

0.3 

0.4 

0.8 

0.3 


3 MATERIAL PROPERTIES AND PARAMETERS USED IN 
THE ANALYSES 


All model predictions that are presented in this work are computed using the 
same set of material data which corresponds to the glass fiber / epoxy material 
used in the experiments described in Section 5. The material data used in the 
analyses is summarized in Table 1. In addition to the elastic properties, intra- 
laminar toughnesses and a slope parameter are required. The toughnesses 
Gi c ,piy and Gn c ,piy are used to compute in-situ strengths (see Section 3.2). 
The slope parameter, p\ 2 = r ] L , determines the slope of the damage activation 
functions for the Mori- Tanaka model and the constant stress model in cr 2 2 — <7i2 
stress space at a 22 = 0, and is chosen as suggested by Puck et al. [26] as 
Pi-2 = Vl = 0.3. 

The damage part of the Mori-Tanaka model requires two additional param- 
eters: the aspect ratio of the Mori-Tanaka inclusions, which is chosen as 
e n = 0.001 to resemble the crack-like geometry, and the damage evolution 
parameter, which is set to k — 0.05 such that damage onset occurs 5% be- 
low the nominal failure load at FI — 0.95. Note that for k — 0, damage 
onset would occur at FI = 1 and the failure index would stay constant during 
damage progression. Hence, for uni-axial transverse tension, the Mori-Tanaka 
model would give the same response as the constant stress model. 


3.1 Plasticity response 


The non-linear shear response of the chosen glass/epoxy material has been 
characterized in Ref. [28] by determining the change of secant shear modulus 
with shear strain based on tensile tests of (±40) s and (±27) s laminates. The 
experimental data shown in Fig. 4b is derived from measurements of laminate’s 
axial modulus, E xx , and Poisson’s ratio, v xy . The results of these tests are used 
to calibrate the plasticity response of the Mori-Tanaka model, the constant 
stress model, and the periodic crack model. The relation between shear stress 
and plastic shear strain (Fig. 4a) can be derived from the original data by 
using the initial shear modulus and factoring out the elastic shear strain. 
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For the periodic crack model, the non-linear shear strain is approximated by 
the function in Eq. (8) with y = 2 x 10 -8 MPa~ 3 [6]. For the Mori-Tanaka 
model, the plastic shear strain is approximated by the analytical curve fit 

with n — 7, k — 147 MPa 

Note that Eq. (13) would be equivalent to the non-linear part of Eq. (8) for 
n — 3. However, it has been found that a higher exponent typically yields a 
better description of the non-linear shear response for composites with epoxy 
or PEEK matrix [10]. 



pi 

712 = 


012 

IT 


The plasticity law for the constant stress model is formulated in terms of the 
normalized secant shear modulus change vs. equivalent shear strain (Fig. 4b). 
Both a linear and a non-linear curve fit are used to define / i2 in Eq. (3) as 

f \ 2 = —28 £ e q + 1-19 . . . linear fit 

/12 = 469 (e eq ) 2 — 38 £eq + 1.22 ... non-linear fit . 

Due to the lack of experimental data for the plastic transverse response, the 
same function is used in the transverse direction, i.e. /2 = fn- Based on the 
non-linear data for a glass fiber/ epoxy material in another study [30], this 
seems to be a reasonable assumption. Nevertheless, this estimation has to be 
kept in mind for predictions that involve plasticity under transverse tension. 




(a) (b) 

Fig. 4. Analytical curve fits for non-linear shear response derived from experimental 
data given in Ref. [28]; (a) Mori-Tanaka and periodic crack model, 072 vs. 7]/, ; (b) 
constant stress model, G^sec/G^ vs. 712- 
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Table 2 

In-situ strengths used in analyses for 4- ply and 8-ply laminates (values in MPa). 


yt y-t 

1 linear 1 non-linear 

8-ply strength 76 72 

4-ply strength 107 93 


S 

73 

82 


3.2 In-situ strength 


Experimental data shows that matrix cracking in an embedded ply occurs 
at higher stresses than in a single ply or UD-laminate (e.g. [31, 32, 33, 34]). 
Assuming that matrix cracks are caused by propagation of pre-existing flaws, 
this ’in-situ’ effect can be explained by the initial flaw size being limited to 
the ply thickness, and the presence of neighboring plies which changes the 
boundary conditions of the underlying fracture mechanics problem. Camanho 
et al. [22] proposed to account for the in-situ effect by introducing in-situ 
strengths that depend on ply thickness. The in-situ strengths are derived based 
on fracture mechanics by comparing the energy released by crack formation 
to the material’s intralaminar fracture toughness. The energy released in the 
embedded ply during crack formation depends on the ply’s constitutive law. 
Therefore, any non-linearity prior to damage onset needs to be accounted for 
in the computation of in-situ strengths. 

Two types of laminates are considered in the current study; laminate A with 
a 12-ply (±/3/904 ) s layup and laminate B with a 19-ply (0/ ± f3< t/0i/2) s layup, 
where the 1/2 subscript indicates a half ply. Damage in these laminates occurs 
primarily in the 8-ply stack of 90°-layers of laminate A and in the angle-plies 
of laminate B with a thickness of four plies for +/3 and — /3, respectively. It 
has been found in a previous study [10], that in-situ strengths need to be used 
for these layers to obtain good correlation with experimental data. However, 
it is not possible to assign individual in-situ strengths for each ply in the cur- 
rently implemented versions of some of the models. Therefore, the same in-situ 
strengths are used for all plies of one laminate in order to produce comparable 
results, i.e. all plies in laminate A are assigned 8-ply in-situ strengths while 
the plies in laminate B have 4-ply in-situ strengths (see Table 2). The effect 
of different in-situ strengths has been studied in detail in [10]. 

For transverse tension, all models except the constant stress model assume 
a linear response and the corresponding in-situ strength V][ near is computed 
as proposed by Camanho et al. [22] from the intralaminar Mode I fracture 
toughness, Gi C)P i y . An in-situ transverse strength for non-linear transverse be- 
havior, 4'n 0n _ linear , is used for the constant stress model following a modification 
of Camanho’s approach by Pinho et al. [12], 
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Table 3 

Non-standard input parameters required by the damage models for predicting trans- 
verse matrix cracking (values indicated by *) are derived from experimental data). 


model 

plasticity 

damage initiation 

damage evolution 

Mori- Tanaka 

n*\k 

Pu 

e n , k 

constant stress 

fl2*\f2 

Vl 

N/A 

periodic crack 


Vl 

N/A 

exp. softening 

N/A 

Vl 

N/A 

ABAQUS 

N/A 

N/A 

N/A 


The in-situ shear strength is computed from the mode II fracture toughness, 
Gii C)P iy, for the non-linear shear formulations of the Mori-Tanaka, the constant 
stress, and the periodic crack models and yields the same value for all three 
models. The exponential softening model and the ABAQUS model use a linear 
shear response. The in-situ shear strength for a linear shear response would 
amount to 101 MPa and 143 MPa for 8-ply and 4-ply stacks, respectively. These 
values are unrealistically high considering the maximum attainable shear stress 
shown in Fig. 4a. Therefore, the in-situ shear strength for non-linear shear is 
used for the last two models as well. Note, however, that this inaccuracy has 
little effect on the predictions presented in the following for the exponential 
softening model and the ABAQUS model since those two models are only applied 
to test cases where damage is dominated by transverse tension. 


3.3 Parameter identification 


Most of the damage models discussed in this work require some input param- 
eters that are not readily available from standard test methods. The param- 
eters required by each of the models in addition to ply stiffnesses, strengths 
and intra-laminar toughnesses are listed in Table 3. Except for the empirical 
plasticity parameters for in-plane shear (marked by *)), all of the parameters 
in Table 3 had to be estimated for the glass fiber material considered here. 
This seems to be a drawback of these models. However, some bounds for an 
appropriate choice of the lacking parameters can be given. 

Plasticity under transverse tension has not yet received much attention in the 
literature. But from a theoretical point of view it is reasonable to set / 2 = fn 
since plasticity under transverse tension is caused by the same mechanism as 
shear plasticity. 

Values for slope parameters, p\ 2 or ? ?l, of typical glass and carbon fiber mate- 
rials can be found in the literature [13, 26] and are consistently in a range of 
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0.25 - 0.35. Note that the slope parameters can also be set to zero, in which 
case the Puck and LaRC 04 failure criteria would yield the same result as the 
simpler Hashin criterion used by the ABAQUS model. 

The Mori-Tanaka model offers two parameters for adjusting the degradation 
behavior. The aspect ratio, e n , has to be chosen very small to resemble the 
crack-like geometry. However, reducing the aspect ratio to values smaller than 
0.001 hardly changes the predicted degradation further [35]. The evolution 
parameter, k, accounts for the statistical variation in the local strength of a 
ply and can therefore be estimated based on the scatter of experimental data. 
Choosing k = 0 corresponds to zero strength variation and predicts the stress 
to remain constant after damage onset, similarly to the constant stress model. 


4 QUALITATIVE COMPARISON OF DAMAGE MODELS 

In this section, the five models are compared qualitatively. The comparison 
emphasizes the differences between the model formulations and implications 
of the modeling assumptions. 


4-1 Plasticity 


Some recent experimental investigations have shown that laminates subjected 
to loads that mainly lead to ply shear stresses exhibit a significant amount of 
residual strains after unloading [36, 37]. Figure 5a shows experimental data of 
a uniaxial tension test on a symmetric ±45 laminate including several loading 
cycles [37] . The loops reveal that the non-linearity caused by the shear stresses 
is initially related to residual strains while the shear modulus (indicated by 



Shear stress-strain curve for cyclic [*45“/-45 u ] ?1 test IH2 



(a) (b) 


Fig. 5. (a) Non-linear shear stress-shear strain response of a symmetric ±45 lam- 
inate under uniaxial tension; (experimental results [37]); (b) schematic of model 
predictions (DO m damage onset). 
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dashed lines) remains approximately constant at first. At higher loads, the 
shear modulus starts to change significantly accompanied by a further accu- 
mulation of residual strain. Note that the hysteresis loops also found in the 
experiments, which could be the result of viscous effects, are not considered 
by any of the models. 

The Mori-Tanaka and the constant stress models account for residual strains 
using plasticity formulations. A schematic of the shear response predicted by 
those two models is shown in Fig. 5b as red solid and green dashed lines, re- 
spectively. Both models assume that the non-linearity prior to damage onset 
(DO) is caused by plasticity only. Therefore, any unloading curves below dam- 
age onset follow the initial slope of G°. After damage onset, the constant stress 
model assumes no further accumulation of plastic strain, i.e. all non-linearity 
after damage onset is caused by matrix damage. In the Mori-Tanaka model, 
the plastic strain is a function of shear stress and since the stress continues to 
increase slightly after damage onset, the amount of plastic strain continues to 
grow in this model. Note that for k — (1, the failure index in the Mori-Tanaka 
model would remain constant during damage progression, which means that 
the stresses and the plastic strain would not increase after damage onset (for 
a constant ratio of shear vs. transverse stress). 

Also shown in Fig. 5b is a schematic of the response from the periodic crack 
model, which accounts for shear non-linearity by a non-linear elastic stress- 
strain law. Therefore, the same non-linear curve is followed for loading and 
unloading as long as the maximum stress does not exceed the damage onset 
stress. After damage onset, the stresses decrease with increasing load and the 
unloading path would be linear and leading back to the origin. 

When only monotonic loading (in each material point) is considered, the defi- 
nition of the unloading path is without consequence. For non-monotonic load- 
ing, which can occur locally due to damage even if the global load is applied 
monotonically, using the correct unloading path is important for predicting 
the redistribution of stresses accurately. 


4-2 Distributed vs. localized cracking 


If a ply embedded in a laminate is subjected to a homogeneous stress or strain 
field, the first matrix crack that spans the whole thickness of a ply typically 
does not lead to laminate failure since stresses can locally be transferred to 
adjacent layers. Consequently, an array of matrix cracks will develop and ac- 
cumulate more or less evenly throughout the ply (left side of Fig. 6). 

The response of a ply containing an array of cracks can be modeled by contin- 
uum mechanics provided that the size of the volume considered is large com- 
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pared to the width of a crack (homogenization requirement). The stresses given 
by a continuum model should equal the volume average of the actual stresses 
(‘stress distribution’ in Fig. 6). As the number of ply cracks increases with 
load, the volume averaged stresses may decrease but as long as the ply inter- 
faces and the adjacent plies are intact, they remain non-zero due to the stresses 
transferred through the ply interfaces (‘stressustrain response’ in Fig. 6). At 
some point, the increase of crack density with load will saturate since the 
stress transfer is limited by the strength of the ply interfaces. The first three 
models presented in Section 2 are intended to model this behavior up to crack 
saturation, with damage progression corresponding to the increasing number 
of ply cracks. Note that delamination after crack saturation can be handled 
by the constant stress model but is not considered in the present study. 

The last two models of Section 2 ( ExpSoft and ABQ) are intended for modeling 
localized damage, i.e. the response due to a single crack. In those models, the 
progression of damage has to be interpreted as the formation of micro-cracks 
and their ultimate coalescence into a ply crack (Fig. 7). For a single ply, it 
is clear that in such a model the stresses have to be reduced to zero in the 
ultimate damage state where the material has been completely separated. For 
a ply embedded in a laminate, ply stresses will be zero at the crack but not 
further away from it as long as the ply interfaces are intact. Consequently, the 
assumption inherent to the localized damage models that shear and transverse 
stresses become zero when the ply crack has fully formed are only valid if ho- 
mogenization is performed over a volume approximately the size of the cracked 
domain (vol . 1 in Fig. 7). The element size in an FEM computation, however, 
will typically be much larger than this for reasons of numerical efficiency (for 
example, similar to vol . 2 in Fig. 7). The stresses averaged over vol . 2 do not 


homogenized volume 


array of cracks 



Fig. 6. Homogenization volume, stresses distribution, and stress-strain response for 
distributed cracking. 
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Fig. 7. Homogenization volumes, stresses distributions, and stress-strain responses 
for localized cracking. 
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tend to zero since a considerable part of the volume is undamaged and loaded 
through shear stresses along the ply interfaces. From these considerations it 
is obvious that the two models for localized damage will not give very accu- 
rate predictions for homogeneous loading since the ply interfaces are expected 
to stay intact until crack saturation. Laminates subjected to inhomogeneous 
stress fields (e.g. near a notch) often develop one dominant ‘splitting’ crack 
accompanied by local delamination (e.g. [38]), in which case stresses vanish in 
a larger area around the crack. 

The present review focuses on damage due to matrix cracking under homoge- 
neous stress / strain states prior to crack saturation from loads that lead to an 
array of matrix cracks but leave the interface intact. As detailed above, the 
localizing damage models are not intended for this purpose. However, they 
will be included in the present comparison in order to point out and discuss 
the effects of different implementation strategies. Analyses using the localized 
damage models ( ExpSoft and ABQ ) are performed by applying the homoge- 
neous load to a single, 4-noded finite element. Since the laminate response 
predicted by the localized damage models is not fixed, but depends on the 
characteristic element length, the element size is chosen as the distance be- 
tween two cracks at crack saturation, 2 L sat , which is estimated by the following 
approach. 


Crack saturation occurs when the distance between two cracks is so small that 
the ply stresses, transmitted via the interface from neighboring plies, do not 
exceed the tensile failure stress. The stresses transmitted through the interface 
are limited by the interface shear strength. Based on a shear lag model for an 
embedded ply under uniaxial transverse tension (Fig. 8a) and assuming the 
shear stress along the ply interface to be constant and equal to the interface 
shear strength, iS IF , the crack spacing at saturation can be estimated as 

2L sat = tY t /S w , (15) 



Fig. 8. (a) Shear lag model used to estimate the crack spacing at crack satura- 
tion, (b) schematic of stress-strain response and dissipated energy density at crack 
saturation. 
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where t denotes the thickness of the cracked ply (or cluster of like oriented 
plies). Based on the linear in-situ strengths for transverse tension (Table 2) 
and estimating the interlaminar shear strength as S ' [ ‘ = 73 MPa, the crack 
spacing is computed from Eq. (15) as 2L sat = 1.2 mm and 2L sat = 0.84 mm for 
the laminates with 8-ply and 4-ply in-situ strengths, respectively. The dissi- 
pated energy density predicted by the localized damage models then amounts 
to U sat = Gi c>p i y /(2L sat ), which equals the energy density dissipated by the 
distributed damage models up to the saturation stage (see Fig. 8b). 

The stress-strain response for an embedded ply under transverse tension is 
shown in Fig. 9a. The experimental results are derived from data provided in 
[29], VarnaOl, and [28], Varna99, for uniaxial tension tests on flat specimens 
with a (02/904) s and a (0/ 90g /0i/ 2 ) s lay-up, respectively, by assuming that the 
change in the laminate response is only the result of a change in the transverse 
Young’s modulus of the 90-ply clusters. Note that a change of Gi 2 is irrelevant 
to the response of a 0/90-laminate under uniaxial tension. 

As can be expected, Fig. 9a shows that the prediction of the periodic crack 
model is the most accurate since that model is an exact solution of the corre- 
sponding boundary value problem. The two localized damage models ExpSoft 
and ABQ predict the stresses to reduce to zero which does not agree with 
the experimental data due to the reasons discussed above. The Mori-Tanaka 
model predicts ply stresses to be almost constant after damage onset, which 
is a result of the damage evolution law chosen. The response of the constant 
stress model is very similar, but with a small amount of non-linearity prior 
to damage onset, due to the plasticity law for transverse tension, and with 
the stress remaining constant after damage onset as imposed by the model’s 
assumption. Note that according to Joffe et al. [27], no cracks were observed 
in the (0 2 /904) s laminate below a strain of e xx 0.6% (e xx — e 22 in Fig. 9). 



Fig. 9. Embedded 90°-ply in 0/90 lay-up under uniaxial tension: (a) stress strain 
response, (b) accumulation of matrix cracks; experimental data from [28, 29]. 


19 


However, the experimental data shows some small non-linearity starting at 
slightly lower strain. This could be an indication that plasticity indeed occurs 
in the 90° ply prior to damage onset, as assumed in the constant stress model. 

Figure 9b shows the increase of crack density, CD, as a function of axial strain 
for the same two cross-ply laminates. In the periodic crack model, crack density 
is an internal variable and is therefore directly available for plotting. Note 
however, that the periodic crack model is strictly only valid for the (02/904) s 
laminate (corresponding to data points labeled exp. VarnaOl [29]) and not to 
the other lay-up. The crack densities for the Mori- Tanaka and constant stress 
models are computed by dividing the predicted energy density dissipated due 
to damage, U dls , with the energy dissipated by the formation of one crack, 
f?ic,piy> he. CD = U dis /Gl C) p\y. The predictions of these two models are the 
same for both laminates (solid red and dashed green lines, respectively, in 
Fig. 9b). 

Also indicated in Fig. 9 is the strain at which the saturation crack density, 
CD sat — 1/(2 L sat ) = 0.83/mm, based on the shear lag model for the (02/904) s 
laminate, is reached. This saturation crack density strain is £ xx « 1.5% for 
the Mori Tanaka and the constant stress models, and e xx « 1.95% for the 
periodic crack model. As noted above, the energy dissipation per unit volume 
predicted by these three models for CD = 0.83/mm is the same as the energy 
density dissipated by the localized damage models for this lay-up. 


4-3 Damage activation and damage evolution 


The damage models considered in this study use stress based failure envelopes 
to predict damage onset. These failure criteria distinguish between matrix and 
fiber failure for both tension and compression, and accordingly, use separate 
functions for each failure mode. For tensile matrix failure in plane stress, 
which is investigated here, only in-plane shear and transverse tensile stresses 
are relevant. The failure envelopes for matrix tension can thus be plotted in 
0 - 12-022 stress space (see Fig. 10) and show only minor differences between 
the failure criteria. The main difference is the slope of the failure envelope 
at 022 = 0. The Hashin criterion assumes zero slope at that point, while the 
Puck 2D and the Mohr-Coulomb based criteria require an additional material 
parameter to determine the slope, and the LaRC 04 criterion determines the 
slope by the ply’s strengths and toughnesses (see Section 3). Note that the 
envelopes in Fig. 10 are plotted for 8-ply in-situ strengths. The in-situ strength 
for transverse tension used by the constant stress model takes plasticity into 
account and is therefore lower than for the other models (cf. Table 2). 


20 



With respect to the damage activation functions, it should be noted that the 
Hashin and Puck 2D equations only contain stress components to the power 
of one. Therefore, as long as the material response is linear, the failure index 
given by Hashin and Puck 2D scales linearly with load, and the failure load is 
given by cr° — cr/FI. The failure index computed using the LaRC04 and the 
Mohr-Coulomb criteria cannot be interpreted this way since the corresponding 
equations contain first and second order terms of cr 2 2 - The advantages of a 
linear failure index, however, are only relevant for linear analyses (e.g. ‘First 
Ply Failure’ analysis). In a damage analysis, the failure index generally does 
not scale with load since the material response is non-linear. 

As long as a ply is undamaged, the application of stress based damage activa- 
tion functions is straight forward. Once a ply has been damaged, the correct 
stresses to use for evaluating the damage activation functions need to be identi- 
fied. Perhaps the most intuitive solution is to use the homogenized ply stresses 
given by the constitutive model as proposed by the Mori-Tanaka model. The 
disadvantage of this approach is that the constitutive equation of the damage 
model needs to be solved iteratively, since the stress state depends on the 
damage state, and vice versa. 

One way to obviate the need for iteration is to introduce an effective stress 
tensor that is independent of the damage state (Eqs. (4) and (7) for models 
ConStress, Per Crack, and Exp Soft). The drawback with this method is that 
it does not account for the changing stress-strain coupling in the longitudinal 
and transverse directions that is caused by damage. Hence, the damage pre- 
diction of secondary failure modes may be inaccurate. Consider, for example, 
a uni-directional laminate under transverse tension that has been damaged 
due to matrix cracking. Because of the damage, the minor Poisson’s ratio, u 2 \ , 



o 2 2 [MPa] 

Fig. 10. Damage activation functions for tensile matrix failure used by the damage 
models: Puck 2D for MoriT , Mohr-Coulomb for ConStress, LaRC04 for PerCrack 
and ExpSoft, and Hashin for ABQ. 
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approaches zero and the contraction in fiber direction vanishes. The effective 
stress is defined as the stress state that would exist if an undamaged ply were 
subjected to the same strain state (i.e. e 22 > 0, £u ~ 0). Since the minor 
Poisson’s ratio of an intact ply is larger than in a damaged ply, this strain 
state results in a non-zero effective stress in the fiber direction, aff > 0- I 11 
the worst case, depending on the material properties and the degradation laws 
used, this effective stress can result in a predicted fiber failure even though 
the ply is only loaded in the transverse direction. 

The ABAQUS damage model also uses an effective stress tensor to predict dam- 
age initiation. However, this tensor is defined by scaling the components of 
the homogenized stress tensor individually based on the corresponding dam- 
age variable, df, d m , or d s , (see Eq. 10) which circumvents the problem above. 
The effective stress tensor used by the ABAQUS damage model depends on 
the damage state (defined by the damage variables). Nevertheless, there is 
no need for iteration because once damage has initiated, the evolution of the 
damage state in this model is a function of the strain state only. The draw- 
back of using strains to control damage evolution, however, is that due to the 
ply coupling constraints in laminates, stresses and strains do not necessarily 
correlate. Damage progression observed experimentally is typically driven by 
stresses rather than strains. For example, uniaxial tension of a 0/90 laminate 
first causes matrix cracking in the 90°-plies, but with increasing load, cracks 
develop also in the 0°-plies. This is a result of the constraint of the Poisson 
effect imposed by the 90°Tayers that leads to transverse tensile stresses in the 
0°-plies. If this loading scenario is analyzed using the ABAQUS damage model, 
the Hashin criterion correctly predicts damage initiation due to matrix tension 
in the 0°-plies. However, there is no damage evolution in those plies because 
the transverse strain, e 22 1 is negative, indicating compression. Consequently, 
the equivalent displacement for matrix tension, which controls damage evolu- 
tion for this failure mode, is zero and the development of cracks in the 0°-plies 
cannot be captured. 


4-4 Stiffness degradation and anisotropy of damage 


The effect of matrix cracks on the homogenized ply stiffnesses varies in dif- 
ferent directions because of the crack geometry. A damage model for matrix 
cracking thus needs to capture this anisotropy of damage by degrading the var- 
ious components of the elasticity tensor differently. For orthotropic materials 
under plane stress, the material response is determined by four independent 
variables E \ , E 2 , G 12, and U\ 2 . Figure 11 shows the change of those four mate- 
rial parameters, normalized by their respective initial values, as predicted by 
the five damage models for a (62/904)5 laminate under uniaxial tension. 
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In the Mori-Tanaka model (Fig. 11a), the different degradation of the prop- 
erties is a result of the chosen aspect ratio of the Mori-Tanaka inclusions. An 
aspect ratio of e n = 1 would lead the degradation of all properties to be equal 
(isotropic damage). To mimic the effect of thin cracks, the aspect ratio has 
to be chosen very small (e.g. e n = 0.001) which results in a barely noticeable 
increase of E\ and u i2 , and a degradation for E 2 that is different than the 
degradation for G 12 ■ Note that all other models assume E\ and v\ 2 to remain 
constant, which leads to a degradation of the minor Poisson’s ratio which is 
equal to that of the transverse Young’s modulus, ^ 21/^21 = ^/E®. The same 
result would be obtained with the Mori-Tanaka model for perfectly thin cracks 
(i.e. e n 0). 
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Fig. 11. Predicted degradation of elastic properties: (a) Mori-Tanaka model; (b) 
constant stress model; (c) periodic crack model; (d) exponential softening model 
and ABAQUS model. 
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In the periodic crack model (Fig. 11c), the degradation of E 2 and G i2 comes 
from the solution of the boundary value problem of a unit cell containing one 
crack. Note that both the Mori- Tanaka and the periodic crack models predict 
a greater degradation in the transverse Young’s modulus than in the shear 
modulus. This observation generally holds true for other material systems as 
well (glass and carbon fiber / epoxy composites). 

For the constant stress model (Fig. lib), the degradation of E 2 is chosen 
such that the ply stress remains constant after damage onset. Therefore, the 
^-degradation of the constant stress model is very similar to that of the Mori- 
Tanaka model, which also predicts an almost constant ply stress, and slightly 
less than that of the periodic crack model, which shows a softening behavior 
(Fig. 9). The degradation of G \ 2 in the constant stress model is assumed to 
equal that of E 2 . Note that for the constant stress model, the transverse tensile 
stresses in the 90°-ply lead to plasticity prior to damage onset, causing the 
linear degradation between e 22 « 0.5 and £22 = 0.6. The curves plotted for the 
constant stress model thus represent the changes of secant moduli rather than 
those of elastic moduli. As a result of the assumption f i2 — f 22 for the given 
material (Sec. 3), the effect of plasticity in the constant stress model here is 
the same for E 2 and G 12 , however, this is not generally the case. 

In the localized damage models, exponential softening and ABAQUS (Fig. lid), 
the degradations of E 2 and G \ 2 are determined by the degradation functions 
for shear and transverse tension and depend on element size. In general, those 
two models predict a much faster degradation of E 2 and G\ 2 than the first 
three models since they assume the stresses to reach zero in the fully damaged 
state (i.e. E 2 — G\ 2 = 0). The degradation of G\ 2 in the ABAQUS model equals 
the degradation of E 2 (solid black line in Fig. lid). In the exponential softening 
model, the degradation of E 2 and G\ 2 are computed separately as functions 
of the respective strength and fracture toughness for transverse tension and 
shear, as well as the characteristic element length. In some cases, this can 
lead to the shear degradation exceeding the transverse degradation. However, 
in the ultimate damage state, the degradation of the two moduli is the same 
since both are degraded to zero. 


5 COMPARISON TO EXPERIMENTAL DATA 


For a quantitative assessment of the five models, predictions are evaluated 
against experimental data from two series of uni-axial tension tests conducted 
by Varna et al. [28, 29]. For each test, the degradation of the laminate’s ax- 
ial modulus and Poisson’s ratio are compared. The parameters for the glass 
fiber / epoxy material used are listed in Table 1. Residual thermal stresses are 
not taken into account in the present study. The effect of curing stresses for 
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the same set of experimental tests has been studied in [10] using the Mori- 
Tanaka damage model and was found to be small compared to the scatter of 
data. 

Since the stress / strain state in the specimens is homogeneous, the loading 
leads to a series of matrix cracks spread uniformly throughout the whole spec- 
imen. The increasing crack density was also recorded during the experiments 
[28, 29]. The localized damage models from Sections 2.4 and 2.5 are not in- 
tended for modeling progressive cracking under homogeneous loading and their 
predictions cannot be expected to correlate well with the experimental data. 
Nevertheless, these two models are included in the comparison as well in order 
to gain an understanding of how big the deviation is. 


5.1 Laminates type A: embedded 90° plies (±/3/904) s [27, 29] 


The laminates type A have a lay-up of (±/?/904) s with four different lay-up 
angles /3 — 0°, 15°, 30°, 40°. The predictions of analyses compared to exper- 
imental data [29] are shown in Figs. 12-15. Damage in these laminates is 
dominated by matrix cracking in the 90° plies and crack density was moni- 
tored in these layers only during the experiments. However, all damage models 
(except for the periodic crack model which assumes cracking to occur only in 
the 90° layer) also predict damage in the angle-plies. The load at which the 
models predict the onset of cracking in the angle-plies depends on the angle /? 
and on the chosen in-situ strengths. Here, 8-ply in-situ strengths (Table 2) are 
used for all plies. The effect on the predicted load response of using different 
in-situ strengths for each ply depending on the ply’s thickness and location is 
addressed in a previous study based on the Mori- Tanaka damage model [10]. 




(b) 


Fig. 12. Results for lay-up (02/904) s ; (a) laminate’s axial modulus and (b) Poisson’s 
ratio, normalized by their initial value; experimental data from [29]. 
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The figures show that, in general, predictions from the three distributed dam- 
age models are fairly similar, especially as long as cracking only occurs in the 
90° plies, and that the correlation with experiments is quite good. The peri- 
odic crack model predicts slightly more degradation of laminate axial stiffness 
than the Mori-Tanaka and the constant stress models (Figs. 12a-15a) and the 
periodic crack predictions correlate very well with the test data for /3 = 0° 
and f3 = 15°. For (3 — 30° and /3 — 40°, damage in the angle plies leads to 
additional degradation of the laminate axial stiffness which is disregarded by 
the periodic crack model. As a result, predictions of the Mori-Tanaka and con- 
stant stress models are slightly better than the periodic crack model for those 
two lay-ups once cracking starts in the angle plies. 
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Fig. 13. Results for lay-up (±15/904) s ; (a) laminate’s axial modulus and (b) Pois- 
son’s ratio, normalized by their initial value; experimental data from [29]. 




Fig. 14. Results for lay-up (±30/904) s ; (a) laminate’s axial modulus and (b) Pois- 
son’s ratio, normalized by their initial value; experimental data from [29]. 
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Similar observations can be made regarding the degradation of laminate Pois- 
son’s ratio (Figs. 12b-15b). All three models predict a similar response as 
long as damage is restricted to the 90° plies with slightly more degradation 
predicted by the periodic crack model. The onset of angle-ply damage can be 
noticed in the Mori-Tanaka curves as a distinctive kink at e xx pc 1.5 ((3 = 15°), 
e xx « 0.9 {(5 — 30°), and e xx « 0.75 (/3 = 40°). The increased degradation of 
laminate Poisson’s ratio due to angle-ply damage seems to agree well with the 
experiments. The constant stress model also predicts damage in the angle plies 
starting at a similar axial strain but the angle-ply damage has the opposite 
effect on the Poisson’s ratio. The degradation of laminate Poisson’s ratio pre- 
dicted by the constant stress model with angle-ply damage is less than without 
angle-ply damage. This is a result of the assumption that the degradations of 
E 2 and G \2 are equal as illustrated in Fig. 16 by means of the (±40/904) s 
lay-up. If the degradation of G\ 2 is set to 60% of the degradation of E 2 (solid 
red line in Fig. 16) instead of equal to the degradation of E 2 as before (dashed 
green line), then the cracking in the angle plies leads to more degradation 
of the laminate’s Poisson’s ratio (Fig. 16b) while the change has hardly any 
effect on the degradation of the laminate’s axial modulus (Fig. 16a). 

As expected, the two localized damage models do not correlate well with 
the experimental data in Figs. 12-15 since they are not applicable to the 
homogeneous loading conditions. The laminate degradation is overestimated 
since E 2 and G \ 2 for the damaged plies are degraded all the way to zero in both 
models. The leveling of the curves from the ABAQUS model signifies complete 
damage (i.e. E 2 = G \ 2 = 0) of one type of layer. For (3 — 0° and (3 = 15° 
(Figs. 12 and 13), the horizontal parts pertain to a fully damaged 90° layer 
and no damage in the other plies. Figure 14 has two horizontal sections; the 
small first section corresponds to a fully damaged 90° layer and the second 



(b) 



Fig. 15. Results for lay-up (±40/904) s ; (a) laminate’s axial modulus and (b) Pois- 
son’s ratio, normalized by their initial value; experimental data from [29]. 
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horizontal section corresponds to all layers being fully damaged. In Fig. 15, 
damage in the ±40 layers starts before the 90° layer has been fully degraded, 
leaving only one horizontal line when all layers have reached their ultimate 
damage state. Note that for the cases where f3 = 30° and (3 = 40°, there 
is a slight increase of Poisson’s ratio with the onset of damage in the angle 
plies. This is also a result of using the same degradation for E 2 and G \ 2 - The 
responses predicted by the exponential softening model are very similar to the 
ABAQUS model predictions except that the progression of the curves is more 
gradual due to the exponential softening law. 


5.2 Laminates type B: off-axis plies (0/ ± f3i/0i/ 2 ) s [%8] 


The second series of tests was performed on five different laminates with a 
general (0/ ± f3i/0i/ 2 ) s lay-up and angles /3 — 90°, 70°, 55°, 40°, 25° [28]. Non- 
linearity in these laminates originates only from the /3-plies which experience 
various stress ratios, cr 22 /(Ji 2 , depending on the angle /?. For shear dominated 
stress ratios, i.e. (3 = 40° and (3 = 25°, no cracking was observed during the ex- 
periments [28] . These two tests can therefore only be captured by models con- 
taining a non-linear shear formulation such as the Mori- Tanaka , the constant 
stress , and the periodic crack models. For the first three tests ((3 — 90°, 70°, 
and 55°), on the other hand, the periodic crack model cannot be used since 
it does not allow for cracking in more than one layer and is therefore not 
applicable to these lay-ups. 




axial strain, e xx [%] axial strain, e xx [%] 

(a) (b) 

Fig. 16. Effect of changing the ratio of degradation between E 2 and G\ 2 on pre- 
dictions of the Constant Stress model by example of the (±40/904) s lay-up; (a) 
laminate’s axial modulus and (b) Poisson’s ratio, normalized by their initial value; 
experimental data from [29]. 
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Comparisons between the analyses’ predictions and experiments are shown in 
Figs. 17-21. Note that the degradation for /3 = 90°, 70°, and 55° is given in 
[28] as a function of crack density, which can be converted to axial strain by 
analytical curve fits for crack density vs. axial strain which are also provided in 
[28] . All experimental data points shown for those three lay-ups correspond to 
crack densities greater than zero. Note that for (3 — 55°, the relation between 
crack density and axial strain varies significantly from specimen to specimen, 
while the stiffness vs. crack density results are consistent. Therefore, two data 
sets are derived using a lower bound and an upper bound curve to fit the crack 
density vs. strain data. This way, the scatter of experimental data is visualized 
in Fig. 19 by the two data sets for lower bound (lb) and upper bound (ub). 




Fig. 17. Results for lay-up (O/SfOs/Oj^sj ( a ) laminate’s axial modulus and (b) Pois- 
son’s ratio, normalized by their initial value; experimental data from [28]. 




Fig. 18. Results for lay-up (0/ ± 704/0 1 / 2 )s; (a) laminate’s axial modulus and (b) 
Poisson’s ratio, normalized by their initial value; experimental data from [28]. 
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The response for f3 — 90° (Fig. 17) is very similar to the first test of the pre- 
vious test series and confirms the earlier results. In the laminate with /3 = 70° 
(Fig. 18), cracking during experiments starts below the strains at which the 
models predict damage onset (the reduction predicted by the constant stress 
model between e xx « 0.5 and e ^ m 0.8 is due to plasticity). For /3 = 55° on the 
other hand, the models unanimously predict damage onset earlier than was 
observed in the experiments even if the scatter of experimental data is consid- 
ered (Fig. 19, with lower bound (lb) and upper bound (ub) of experimental 
data). The reason for this is not quite clear but could be related to matrix 
plasticity. The test data for the 55° case suggests that there is a significant 
amount of plasticity involved since the axial stiffness has already decreased 
prior to damage onset (first point of each experimental data set in Fig. 19). 




Fig. 19. Results for lay-up (0/ ± 554/0 1 / 2 ) s ; (a) laminate’s axial modulus and (b) 
Poisson’s ratio, normalized by their initial value; experimental data from [28]. 




Fig. 20. Results for lay-up (0/ ± 404/0 1 / 2 ) a ; (a) laminate’s axial modulus and (b) 
Poisson’s ratio, normalized by their initial value; experimental data from [28]. 
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Also the Mori- Tanaka and the constant stress models predict plastic strains 
for that laminate. The interaction between plasticity and damage onset will 
be discussed in more detail in Section 6. 

The predicted change in Poisson’s ratio for (3 — 70° and (3 = 55° is very 
sensitive to the exact model formulations. This sensitivity is best explained 
by reference to the 55° case (Fig. 19b). The Mori-Tanaka model predicts a 
Poisson’s ratio increase due to plasticity, but a decrease due to damage. This 
prediction correlates very well qualitatively with the test data, which also 
shows an increase prior to damage onset (first data point) and slightly de- 
creases thereafter. The plasticity part of the constant stress model predicts a 
very similar response. However, the damage formulation of that model leads 
to a further increase of the Poisson’s ratio, which results from the assumption 
that the degradations of E 2 and G \ 2 are equal. Since the ABAQUS model also 
uses this assumption, it also predicts the Poisson’s ratio to increase with dam- 
age but the effect is much more pronounced since degradation is higher in the 
ABAQUS model. In the exponential softening model, G\ 2 is degraded faster than 
E 2 for the given material and element size, which leads first to a reduction of 
u xy . Ultimately, the ply stiffnesses E 2 and G \ 2 are both reduced to zero, and 
the prediction converges to that of the ABAQUS model. 

The angle plies of the last two lay-ups (/3 = 40° and (3 = 25° in Figs. 20 and 
21, respectively) experience a very small stress ratio a 22 /ai 2 and do not show 
any damage in the strain range considered. The non-linear elastic function 
used by the periodic crack model for the non-linear shear response predicts 
too much non-linearity at low strains, but gives good results at higher strains. 
The predictions of the Mori- Tanaka and the constant sti'ess models are very 
similar up to e xx 1%. Up to that strain, the non-linear shear response is 
known from experiments (see Fig. 4). For s xx > 1%, the plasticity law has to 




(b) 


Fig. 21. Results for lay-up (0/ ± 2bi/t)i/ 2 ) s \ (a) laminate’s axial modulus and (b) 
Poisson’s ratio, normalized by their initial value; experimental data from [28]. 
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be extrapolated. The extrapolation used for the Mori-Tanaka model correlates 
quite well with the experimental data. For the constant stress model, the linear 
extrapolation gives good results for the laminate’s axial stiffness but leads to 
an overprediction of Poisson’s ratio at e xx > 1.5%. Changing to a non-linear 
extrapolation improves the Poisson’s ratio predictions only slightly. 


6 DISCUSSION AND CONCLUSIONS 


In the present study, five different damage models are evaluated for matrix 
dominated loading conditions and are compared both qualitatively and quan- 
titatively. The emphasis is laid on predicting the effects of distributed damage 
in an embedded ply typically caused by fairly homogeneous stress states. Three 
of the damage models considered are specifically designed for damage under 
homogeneous loading, while the other two are localized damage models. Al- 
though the localized damage models are not intended for modeling distributed 
matrix cracking, they are also included in this study for comparison. 

In general, the three distributed damage models correlate well with the ex- 
perimental data shown here. However, none of the models captures all of the 
non-linear effects. Based on comparisons of the models’ predictions to experi- 
mental data of uniaxial tension specimens with various lay-ups, several lessons 
can be learned. 


Non-linearity without damage (plasticity) 

Experimental data shows that the response of a laminate can be non-linear 
prior to matrix cracking, especially under shear dominated loading condi- 
tions. The origin of this non-linear behavior is still disputed but could be due 
to matrix plasticity. The Mori- Tanaka and the constant stress models assume 
plasticity laws while the periodic crack model uses a non-linear elastic formula- 
tion. However, since only monotonic loading was considered here, the validity 
of either plastic or elastic unloading was not investigated. 

All three models for distributed damage capture the non-linear shear response 
based on curve fits of experimental data. Consequently, the agreement be- 
tween predictions and experiments only depends on the quality of the curve 
fit. The expressions used by the Mori-Tanaka and the constant stress mod- 
els capture the non-linear shear response very well, while the periodic crack 
model overpredicts the amount of non-linearity at low strains. Note that the 
correlation of the periodic crack model could easily be improved by using a 
higher exponent in the curve-fit equation. 
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For some of the test cases considered, the non-linear shear response needs to 
be extrapolated. The equations used by the Mori-Tanaka and the periodic 
crack models yield good results, even at high shear strains where extrapo- 
lation is necessary. For the constant stress model, a linear relation between 
secant shear modulus and equivalent strain has been suggested. Although the 
linear approximation matches experimental data at low strains, a non-linear 
extrapolation is needed at higher strains. 

While the Mori-Tanaka and the periodic crack models assume non-linearity 
prior to damage to occur only due to shear stresses, the constant stress model 
also postulates a non-linear response under transverse tension, even though 
non-linearity is not observed in transverse tension tests of UD-laminates. This 
assumption is reasonable for the following reasons. If the non-linearity ob- 
served under shear is a result of plasticity in the matrix, transverse tension 
should also result in plasticity due to the deviatoric part of the strain tensor. 
However, this effect is difficult to measure experimentally, since it only occurs 
in embedded plies as a result of the increased in-situ strength. A UD-laminate 
typically fails at lower stresses than the same ply embedded in a laminate, 
and therefore, does not exhibit any non-linearity prior to failure. Note that 
the test data for laminates with embedded 90° plies show a small amount of 
non-linearity prior to damage onset, which could be an indication of plasticity. 
However, the effect is very small and further study is necessary to draw any 
firm conclusions. 


Damage initiation 

The damage initiation criteria used by the five models for matrix tension are 
very similar. The main difference lies in the definition of the stresses that are 
used to evaluate the initiation functions. The Mori-Tanaka model uses the 
nominal ply stress state, which is a function of the damage state, and requires 
iteration inside the material subroutine to determine the state variables. All 
other models use an effective stress and do not require iteration. However, in 
specific cases this approach can lead to results that are inconsistent with ex- 
perimental observations (depending on the definition of the effective stresses). 

As indicated by the test cases with off-axis plies, the prediction of damage 
onset under multi-axial stress states remains the primary challenge, since none 
of the initiation criteria predict damage onset in the off-axis plies correctly. 
Strains at damage onset are overpredicted for (3 = 70° and underpredicted for 
fd — 55°. The reason for this discrepancy is not clear. A possible explanation 
is the effect of matrix plasticity. When plasticity in the matrix is an issue, the 
stress-strain curve becomes very flat such that a large increase in strain leads 
only to a small increase in stress. As a result, damage initiation is relatively 
insensitive to variations in strain since damage initiation is triggered by stress 
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based criteria. This has two implications. First, it explains the large scatter of 
experimental failure strains for the (5 = 55° test. Second, both the stress state 
and the ply strength need to be known very accurately in order to correctly 
predict damage onset in a model and this requires an adequate plasticity 
formulation and reliable failure criteria for multi-axial stress states. 

Another issue is the effect of plasticity on the perceived strength of an embed- 
ded ply. When the in-situ effect is studied in experiments, only the strain state 
at failure is known and the ply’s failure stress has to be computed based on an 
assumed ply constitutive law. As discussed above, there are several indications 
that an embedded ply experiences plastic deformations, even under transverse 
tensile loading. However, the in-situ strength under transverse tension is typi- 
cally determined based on a linear-elastic transverse response, thereby leading 
to higher transverse stresses and a higher perceived in-situ strength. From 
these considerations, it is clear that the constitutive behavior of a ply has to 
be fully understood, both from the modeling and the experimental perspective, 
in order to improve the prediction of damage onset. 


Strain softening 

The test data of uni-axial tension tests on cross-ply laminates show that dam- 
age under transverse tension leads to strain softening but the ply stresses do 
not vanish completely in the strain range considered. This response is cap- 
tured best by the periodic crack model as long as cracking occurs in one layer 
only (as assumed in the formulation of the model). The Mori-Tanaka and the 
constant stress models assume no softening due to damage which does not 
agree that well with experiments. However, the error in the predicted overall 
response of a multi-axial laminate is small enough to be acceptable for most 
practical purposes. Since a non-softening constitutive law is easier to handle 
from a numerical point of view, the better accuracy of a model with strain 
softening should be weighted against the higher computational effort. 


Property degradation 

The longitudinal Young’s modulus and the major Poisson’s ratio of a ply are 
not affected by matrix cracking if the fibers are perfectly aligned and the 
matrix cracks are perfectly thin. All models considered here generally concur 
in this assumption. In the Mori-Tanaka model, the property degradation is 
controlled by the aspect ratio of the Mori-Tanaka inclusions. If the aspect ratio 
approaches zero, E\ and zq 2 remain constant, independently of the damage 
state. For a slightly higher aspect ratio, E\ is degraded and zq 2 increases with 
damage, which qualitatively may give better correlation with experimental 
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data if the fibers in a ply are not perfectly aligned. However, variation of the 
aspect ratio was not investigated in the present study. 

Based on comparison of model predictions to experimental data of tensile tests 
on various laminates, there is little difference between the predictions of the 
three distributed damage models for the degradation of laminate axial mod- 
ulus. The predicted changes of the laminates’ Poisson’s ratios with damage, 
however, point out some differences in the model formulations. The predicted 
values of Poisson’s ratio qualitatively show better correlation with experiments 
for those models where the degradation of a ply’s transverse Young’s modulus 
is larger than the degradation of the shear modulus ( Mori-Tanaka and peri- 
odic crack). The constant stress model assumes equal degradation of Young’s 
modulus and shear modulus, which can have an effect on the laminate Pois- 
son’s ratio that is opposite to that shown by test data (results for /3 = 30° 
and f3 — 40° in Section 5.1 and fd = 55° in Section 5.2). 

The observation that the Poisson’s ratio is much more sensitive to the different 
model formulations also emphasizes the importance of gathering as much in- 
formation as possible during experimental testing. Looking only at the change 
in axial modulus, there is little difference between predictions from the three 
distributed damage models. The comparisons of Poisson’s ratio, on the other 
hand, give some important clues as to which model assumptions capture the 
actual material response more realistically. 


Localized damage 

The test cases considered in this study are concerned with homogeneous stress 
states that lead to distributed damage. However, failure in composite struc- 
tures is often triggered near stress concentrations which can lead to localized 
cracks, usually combined with delamination near the intersection between the 
crack and the ply interface. The exponential softening model and the ABAQUS 
model are designed to model localized cracking. Since most structures have 
homogeneously loaded areas and geometric features resulting in stress concen- 
trations, the optimal damage model would be a constitutive law that is able to 
distinguish between these two situations and predict the response accordingly. 
However, that would require a non-local formulation of the constitutive law 
which is difficult to realize within the framework of commercial FEM software. 

One alternative is to use models for localized damage throughout the struc- 
ture even though their predictions in homogeneously loaded areas will be less 
accurate. Applying the two localized damage models to the test cases in this 
study shows that the effect of matrix damage in general is overpredicted and, 
depending on the laminate lay-up, the error can be quite high. In addition, 
the exponential softening model and the ABAQUS model do not include any 
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plasticity formulations and therefore do not capture the response under shear 
dominated loading very well. However, the laminates studied here were de- 
signed to emphasize the effect of matrix damage such that it could be studied 
properly. Laminates for practical applications usually have more evenly dis- 
tributed fiber directions (e.g. quasi-isotropic laminates) and thinner ply groups 
so the error in the predicted laminate stiffness would be smaller. Furthermore, 
the predictions can be considered a conservative estimate. It should also be 
kept in mind that the main purpose of these models is to predict ultimate fail- 
ure. In practical applications, especially when quasi-isotropic laminates are 
used, ultimate failure is typically dominated by fiber failure while the exact 
matrix response is often less relevant. 
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